home *** CD-ROM | disk | FTP | other *** search
/ Cream of the Crop 26 / Cream of the Crop 26.iso / os2 / octa209s.zip / octave-2.09 / libcruft / dassl / ddassl.f < prev    next >
Text File  |  1996-07-19  |  65KB  |  1,605 lines
  1.       SUBROUTINE DDASSL (RES, NEQ, T, Y, YPRIME, TOUT, INFO, RTOL, ATOL,
  2.      +   IDID, RWORK, LRW, IWORK, LIW, RPAR, IPAR, JAC)
  3. C***BEGIN PROLOGUE  DDASSL
  4. C***PURPOSE  This code solves a system of differential/algebraic
  5. C            equations of the form G(T,Y,YPRIME) = 0.
  6. C***LIBRARY   SLATEC (DASSL)
  7. C***CATEGORY  I1A2
  8. C***TYPE      DOUBLE PRECISION (SDASSL-S, DDASSL-D)
  9. C***KEYWORDS  DIFFERENTIAL/ALGEBRAIC, BACKWARD DIFFERENTIATION FORMULAS,
  10. C             IMPLICIT DIFFERENTIAL SYSTEMS
  11. C***AUTHOR  PETZOLD, LINDA R., (LLNL)
  12. C             COMPUTING AND MATHEMATICS RESEARCH DIVISION
  13. C             LAWRENCE LIVERMORE NATIONAL LABORATORY
  14. C             L - 316, P.O. BOX 808,
  15. C             LIVERMORE, CA.    94550
  16. C***DESCRIPTION
  17. C
  18. C *Usage:
  19. C
  20. C      EXTERNAL RES, JAC
  21. C      INTEGER NEQ, INFO(N), IDID, LRW, LIW, IWORK(LIW), IPAR
  22. C      DOUBLE PRECISION T, Y(NEQ), YPRIME(NEQ), TOUT, RTOL, ATOL,
  23. C     *   RWORK(LRW), RPAR
  24. C
  25. C      CALL DDASSL (RES, NEQ, T, Y, YPRIME, TOUT, INFO, RTOL, ATOL,
  26. C     *   IDID, RWORK, LRW, IWORK, LIW, RPAR, IPAR, JAC)
  27. C
  28. C
  29. C *Arguments:
  30. C  (In the following, all real arrays should be type DOUBLE PRECISION.)
  31. C
  32. C  RES:EXT     This is a subroutine which you provide to define the
  33. C              differential/algebraic system.
  34. C
  35. C  NEQ:IN      This is the number of equations to be solved.
  36. C
  37. C  T:INOUT     This is the current value of the independent variable.
  38. C
  39. C  Y(*):INOUT  This array contains the solution components at T.
  40. C
  41. C  YPRIME(*):INOUT  This array contains the derivatives of the solution
  42. C              components at T.
  43. C
  44. C  TOUT:IN     This is a point at which a solution is desired.
  45. C
  46. C  INFO(N):IN  The basic task of the code is to solve the system from T
  47. C              to TOUT and return an answer at TOUT.  INFO is an integer
  48. C              array which is used to communicate exactly how you want
  49. C              this task to be carried out.  (See below for details.)
  50. C              N must be greater than or equal to 15.
  51. C
  52. C  RTOL,ATOL:INOUT  These quantities represent relative and absolute
  53. C              error tolerances which you provide to indicate how
  54. C              accurately you wish the solution to be computed.  You
  55. C              may choose them to be both scalars or else both vectors.
  56. C              Caution:  In Fortran 77, a scalar is not the same as an
  57. C                        array of length 1.  Some compilers may object
  58. C                        to using scalars for RTOL,ATOL.
  59. C
  60. C  IDID:OUT    This scalar quantity is an indicator reporting what the
  61. C              code did.  You must monitor this integer variable to
  62. C              decide  what action to take next.
  63. C
  64. C  RWORK:WORK  A real work array of length LRW which provides the
  65. C              code with needed storage space.
  66. C
  67. C  LRW:IN      The length of RWORK.  (See below for required length.)
  68. C
  69. C  IWORK:WORK  An integer work array of length LIW which probides the
  70. C              code with needed storage space.
  71. C
  72. C  LIW:IN      The length of IWORK.  (See below for required length.)
  73. C
  74. C  RPAR,IPAR:IN  These are real and integer parameter arrays which
  75. C              you can use for communication between your calling
  76. C              program and the RES subroutine (and the JAC subroutine)
  77. C
  78. C  JAC:EXT     This is the name of a subroutine which you may choose
  79. C              to provide for defining a matrix of partial derivatives
  80. C              described below.
  81. C
  82. C  Quantities which may be altered by DDASSL are:
  83. C     T, Y(*), YPRIME(*), INFO(1), RTOL, ATOL,
  84. C     IDID, RWORK(*) AND IWORK(*)
  85. C
  86. C *Description
  87. C
  88. C  Subroutine DDASSL uses the backward differentiation formulas of
  89. C  orders one through five to solve a system of the above form for Y and
  90. C  YPRIME.  Values for Y and YPRIME at the initial time must be given as
  91. C  input.  These values must be consistent, (that is, if T,Y,YPRIME are
  92. C  the given initial values, they must satisfy G(T,Y,YPRIME) = 0.).  The
  93. C  subroutine solves the system from T to TOUT.  It is easy to continue
  94. C  the solution to get results at additional TOUT.  This is the interval
  95. C  mode of operation.  Intermediate results can also be obtained easily
  96. C  by using the intermediate-output capability.
  97. C
  98. C  The following detailed description is divided into subsections:
  99. C    1. Input required for the first call to DDASSL.
  100. C    2. Output after any return from DDASSL.
  101. C    3. What to do to continue the integration.
  102. C    4. Error messages.
  103. C
  104. C
  105. C  -------- INPUT -- WHAT TO DO ON THE FIRST CALL TO DDASSL ------------
  106. C
  107. C  The first call of the code is defined to be the start of each new
  108. C  problem. Read through the descriptions of all the following items,
  109. C  provide sufficient storage space for designated arrays, set
  110. C  appropriate variables for the initialization of the problem, and
  111. C  give information about how you want the problem to be solved.
  112. C
  113. C
  114. C  RES -- Provide a subroutine of the form
  115. C             SUBROUTINE RES(T,Y,YPRIME,DELTA,IRES,RPAR,IPAR)
  116. C         to define the system of differential/algebraic
  117. C         equations which is to be solved. For the given values
  118. C         of T,Y and YPRIME, the subroutine should
  119. C         return the residual of the defferential/algebraic
  120. C         system
  121. C             DELTA = G(T,Y,YPRIME)
  122. C         (DELTA(*) is a vector of length NEQ which is
  123. C         output for RES.)
  124. C
  125. C         Subroutine RES must not alter T,Y or YPRIME.
  126. C         You must declare the name RES in an external
  127. C         statement in your program that calls DDASSL.
  128. C         You must dimension Y,YPRIME and DELTA in RES.
  129. C
  130. C         IRES is an integer flag which is always equal to
  131. C         zero on input. Subroutine RES should alter IRES
  132. C         only if it encounters an illegal value of Y or
  133. C         a stop condition. Set IRES = -1 if an input value
  134. C         is illegal, and DDASSL will try to solve the problem
  135. C         without getting IRES = -1. If IRES = -2, DDASSL
  136. C         will return control to the calling program
  137. C         with IDID = -11.
  138. C
  139. C         RPAR and IPAR are real and integer parameter arrays which
  140. C         you can use for communication between your calling program
  141. C         and subroutine RES. They are not altered by DDASSL. If you
  142. C         do not need RPAR or IPAR, ignore these parameters by treat-
  143. C         ing them as dummy arguments. If you do choose to use them,
  144. C         dimension them in your calling program and in RES as arrays
  145. C         of appropriate length.
  146. C
  147. C  NEQ -- Set it to the number of differential equations.
  148. C         (NEQ .GE. 1)
  149. C
  150. C  T -- Set it to the initial point of the integration.
  151. C         T must be defined as a variable.
  152. C
  153. C  Y(*) -- Set this vector to the initial values of the NEQ solution
  154. C         components at the initial point. You must dimension Y of
  155. C         length at least NEQ in your calling program.
  156. C
  157. C  YPRIME(*) -- Set this vector to the initial values of the NEQ
  158. C         first derivatives of the solution components at the initial
  159. C         point.  You must dimension YPRIME at least NEQ in your
  160. C         calling program. If you do not know initial values of some
  161. C         of the solution components, see the explanation of INFO(11).
  162. C
  163. C  TOUT -- Set it to the first point at which a solution
  164. C         is desired. You can not take TOUT = T.
  165. C         integration either forward in T (TOUT .GT. T) or
  166. C         backward in T (TOUT .LT. T) is permitted.
  167. C
  168. C         The code advances the solution from T to TOUT using
  169. C         step sizes which are automatically selected so as to
  170. C         achieve the desired accuracy. If you wish, the code will
  171. C         return with the solution and its derivative at
  172. C         intermediate steps (intermediate-output mode) so that
  173. C         you can monitor them, but you still must provide TOUT in
  174. C         accord with the basic aim of the code.
  175. C
  176. C         The first step taken by the code is a critical one
  177. C         because it must reflect how fast the solution changes near
  178. C         the initial point. The code automatically selects an
  179. C         initial step size which is practically always suitable for
  180. C         the problem. By using the fact that the code will not step
  181. C         past TOUT in the first step, you could, if necessary,
  182. C         restrict the length of the initial step size.
  183. C
  184. C         For some problems it may not be permissible to integrate
  185. C         past a point TSTOP because a discontinuity occurs there
  186. C         or the solution or its derivative is not defined beyond
  187. C         TSTOP. When you have declared a TSTOP point (SEE INFO(4)
  188. C         and RWORK(1)), you have told the code not to integrate
  189. C         past TSTOP. In this case any TOUT beyond TSTOP is invalid
  190. C         input.
  191. C
  192. C  INFO(*) -- Use the INFO array to give the code more details about
  193. C         how you want your problem solved.  This array should be
  194. C         dimensioned of length 15, though DDASSL uses only the first
  195. C         eleven entries.  You must respond to all of the following
  196. C         items, which are arranged as questions.  The simplest use
  197. C         of the code corresponds to answering all questions as yes,
  198. C         i.e. setting all entries of INFO to 0.
  199. C
  200. C       INFO(1) - This parameter enables the code to initialize
  201. C              itself. You must set it to indicate the start of every
  202. C              new problem.
  203. C
  204. C          **** Is this the first call for this problem ...
  205. C                Yes - Set INFO(1) = 0
  206. C                 No - Not applicable here.
  207. C                      See below for continuation calls.  ****
  208. C
  209. C       INFO(2) - How much accuracy you want of your solution
  210. C              is specified by the error tolerances RTOL and ATOL.
  211. C              The simplest use is to take them both to be scalars.
  212. C              To obtain more flexibility, they can both be vectors.
  213. C              The code must be told your choice.
  214. C
  215. C          **** Are both error tolerances RTOL, ATOL scalars ...
  216. C                Yes - Set INFO(2) = 0
  217. C                      and input scalars for both RTOL and ATOL
  218. C                 No - Set INFO(2) = 1
  219. C                      and input arrays for both RTOL and ATOL ****
  220. C
  221. C       INFO(3) - The code integrates from T in the direction
  222. C              of TOUT by steps. If you wish, it will return the
  223. C              computed solution and derivative at the next
  224. C              intermediate step (the intermediate-output mode) or
  225. C              TOUT, whichever comes first. This is a good way to
  226. C              proceed if you want to see the behavior of the solution.
  227. C              If you must have solutions at a great many specific
  228. C              TOUT points, this code will compute them efficiently.
  229. C
  230. C          **** Do you want the solution only at
  231. C                TOUT (and not at the next intermediate step) ...
  232. C                 Yes - Set INFO(3) = 0
  233. C                  No - Set INFO(3) = 1 ****
  234. C
  235. C       INFO(4) - To handle solutions at a great many specific
  236. C              values TOUT efficiently, this code may integrate past
  237. C              TOUT and interpolate to obtain the result at TOUT.
  238. C              Sometimes it is not possible to integrate beyond some
  239. C              point TSTOP because the equation changes there or it is
  240. C              not defined past TSTOP. Then you must tell the code
  241. C              not to go past.
  242. C
  243. C           **** Can the integration be carried out without any
  244. C                restrictions on the independent variable T ...
  245. C                 Yes - Set INFO(4)=0
  246. C                  No - Set INFO(4)=1
  247. C                       and define the stopping point TSTOP by
  248. C                       setting RWORK(1)=TSTOP ****
  249. C
  250. C       INFO(5) - To solve differential/algebraic problems it is
  251. C              necessary to use a matrix of partial derivatives of the
  252. C              system of differential equations. If you do not
  253. C              provide a subroutine to evaluate it analytically (see
  254. C              description of the item JAC in the call list), it will
  255. C              be approximated by numerical differencing in this code.
  256. C              although it is less trouble for you to have the code
  257. C              compute partial derivatives by numerical differencing,
  258. C              the solution will be more reliable if you provide the
  259. C              derivatives via JAC. Sometimes numerical differencing
  260. C              is cheaper than evaluating derivatives in JAC and
  261. C              sometimes it is not - this depends on your problem.
  262. C
  263. C           **** Do you want the code to evaluate the partial
  264. C                derivatives automatically by numerical differences ...
  265. C                   Yes - Set INFO(5)=0
  266. C                    No - Set INFO(5)=1
  267. C                  and provide subroutine JAC for evaluating the
  268. C                  matrix of partial derivatives ****
  269. C
  270. C       INFO(6) - DDASSL will perform much better if the matrix of
  271. C              partial derivatives, DG/DY + CJ*DG/DYPRIME,
  272. C              (here CJ is a scalar determined by DDASSL)
  273. C              is banded and the code is told this. In this
  274. C              case, the storage needed will be greatly reduced,
  275. C              numerical differencing will be performed much cheaper,
  276. C              and a number of important algorithms will execute much
  277. C              faster. The differential equation is said to have
  278. C              half-bandwidths ML (lower) and MU (upper) if equation i
  279. C              involves only unknowns Y(J) with
  280. C                             I-ML .LE. J .LE. I+MU
  281. C              for all I=1,2,...,NEQ. Thus, ML and MU are the widths
  282. C              of the lower and upper parts of the band, respectively,
  283. C              with the main diagonal being excluded. If you do not
  284. C              indicate that the equation has a banded matrix of partial
  285. C              derivatives, the code works with a full matrix of NEQ**2
  286. C              elements (stored in the conventional way). Computations
  287. C              with banded matrices cost less time and storage than with
  288. C              full matrices if 2*ML+MU .LT. NEQ. If you tell the
  289. C              code that the matrix of partial derivatives has a banded
  290. C              structure and you want to provide subroutine JAC to
  291. C              compute the partial derivatives, then you must be careful
  292. C              to store the elements of the matrix in the special form
  293. C              indicated in the description of JAC.
  294. C
  295. C          **** Do you want to solve the problem using a full
  296. C               (dense) matrix (and not a special banded
  297. C               structure) ...
  298. C                Yes - Set INFO(6)=0
  299. C                 No - Set INFO(6)=1
  300. C                       and provide the lower (ML) and upper (MU)
  301. C                       bandwidths by setting
  302. C                       IWORK(1)=ML
  303. C                       IWORK(2)=MU ****
  304. C
  305. C
  306. C        INFO(7) -- You can specify a maximum (absolute value of)
  307. C              stepsize, so that the code
  308. C              will avoid passing over very
  309. C              large regions.
  310. C
  311. C          ****  Do you want the code to decide
  312. C                on its own maximum stepsize?
  313. C                Yes - Set INFO(7)=0
  314. C                 No - Set INFO(7)=1
  315. C                      and define HMAX by setting
  316. C                      RWORK(2)=HMAX ****
  317. C
  318. C        INFO(8) -- Differential/algebraic problems
  319. C              may occaisionally suffer from
  320. C              severe scaling difficulties on the
  321. C              first step. If you know a great deal
  322. C              about the scaling of your problem, you can
  323. C              help to alleviate this problem by
  324. C              specifying an initial stepsize HO.
  325. C
  326. C          ****  Do you want the code to define
  327. C                its own initial stepsize?
  328. C                Yes - Set INFO(8)=0
  329. C                 No - Set INFO(8)=1
  330. C                      and define HO by setting
  331. C                      RWORK(3)=HO ****
  332. C
  333. C        INFO(9) -- If storage is a severe problem,
  334. C              you can save some locations by
  335. C              restricting the maximum order MAXORD.
  336. C              the default value is 5. for each
  337. C              order decrease below 5, the code
  338. C              requires NEQ fewer locations, however
  339. C              it is likely to be slower. In any
  340. C              case, you must have 1 .LE. MAXORD .LE. 5
  341. C          ****  Do you want the maximum order to
  342. C                default to 5?
  343. C                Yes - Set INFO(9)=0
  344. C                 No - Set INFO(9)=1
  345. C                      and define MAXORD by setting
  346. C                      IWORK(3)=MAXORD ****
  347. C
  348. C        INFO(10) --If you know that the solutions to your equations
  349. C               will always be nonnegative, it may help to set this
  350. C               parameter. However, it is probably best to
  351. C               try the code without using this option first,
  352. C               and only to use this option if that doesn't
  353. C               work very well.
  354. C           ****  Do you want the code to solve the problem without
  355. C                 invoking any special nonnegativity constraints?
  356. C                  Yes - Set INFO(10)=0
  357. C                   No - Set INFO(10)=1
  358. C
  359. C        INFO(11) --DDASSL normally requires the initial T,
  360. C               Y, and YPRIME to be consistent. That is,
  361. C               you must have G(T,Y,YPRIME) = 0 at the initial
  362. C               time. If you do not know the initial
  363. C               derivative precisely, you can let DDASSL try
  364. C               to compute it.
  365. C          ****   Are the initialHE INITIAL T, Y, YPRIME consistent?
  366. C                 Yes - Set INFO(11) = 0
  367. C                  No - Set INFO(11) = 1,
  368. C                       and set YPRIME to an initial approximation
  369. C                       to YPRIME.  (If you have no idea what
  370. C                       YPRIME should be, set it to zero. Note
  371. C                       that the initial Y should be such
  372. C                       that there must exist a YPRIME so that
  373. C                       G(T,Y,YPRIME) = 0.)
  374. C
  375. C  RTOL, ATOL -- You must assign relative (RTOL) and absolute (ATOL
  376. C         error tolerances to tell the code how accurately you
  377. C         want the solution to be computed.  They must be defined
  378. C         as variables because the code may change them.  You
  379. C         have two choices --
  380. C               Both RTOL and ATOL are scalars. (INFO(2)=0)
  381. C               Both RTOL and ATOL are vectors. (INFO(2)=1)
  382. C         in either case all components must be non-negative.
  383. C
  384. C         The tolerances are used by the code in a local error
  385. C         test at each step which requires roughly that
  386. C               ABS(LOCAL ERROR) .LE. RTOL*ABS(Y)+ATOL
  387. C         for each vector component.
  388. C         (More specifically, a root-mean-square norm is used to
  389. C         measure the size of vectors, and the error test uses the
  390. C         magnitude of the solution at the beginning of the step.)
  391. C
  392. C         The true (global) error is the difference between the
  393. C         true solution of the initial value problem and the
  394. C         computed approximation.  Practically all present day
  395. C         codes, including this one, control the local error at
  396. C         each step and do not even attempt to control the global
  397. C         error directly.
  398. C         Usually, but not always, the true accuracy of the
  399. C         computed Y is comparable to the error tolerances. This
  400. C         code will usually, but not always, deliver a more
  401. C         accurate solution if you reduce the tolerances and
  402. C         integrate again.  By comparing two such solutions you
  403. C         can get a fairly reliable idea of the true error in the
  404. C         solution at the bigger tolerances.
  405. C
  406. C         Setting ATOL=0. results in a pure relative error test on
  407. C         that component.  Setting RTOL=0. results in a pure
  408. C         absolute error test on that component.  A mixed test
  409. C         with non-zero RTOL and ATOL corresponds roughly to a
  410. C         relative error test when the solution component is much
  411. C         bigger than ATOL and to an absolute error test when the
  412. C         solution component is smaller than the threshhold ATOL.
  413. C
  414. C         The code will not attempt to compute a solution at an
  415. C         accuracy unreasonable for the machine being used.  It will
  416. C         advise you if you ask for too much accuracy and inform
  417. C         you as to the maximum accuracy it believes possible.
  418. C
  419. C  RWORK(*) --  Dimension this real work array of length LRW in your
  420. C         calling program.
  421. C
  422. C  LRW -- Set it to the declared length of the RWORK array.
  423. C               You must have
  424. C                    LRW .GE. 40+(MAXORD+4)*NEQ+NEQ**2
  425. C               for the full (dense) JACOBIAN case (when INFO(6)=0), or
  426. C                    LRW .GE. 40+(MAXORD+4)*NEQ+(2*ML+MU+1)*NEQ
  427. C               for the banded user-defined JACOBIAN case
  428. C               (when INFO(5)=1 and INFO(6)=1), or
  429. C                     LRW .GE. 40+(MAXORD+4)*NEQ+(2*ML+MU+1)*NEQ
  430. C                           +2*(NEQ/(ML+MU+1)+1)
  431. C               for the banded finite-difference-generated JACOBIAN case
  432. C               (when INFO(5)=0 and INFO(6)=1)
  433. C
  434. C  IWORK(*) --  Dimension this integer work array of length LIW in
  435. C         your calling program.
  436. C
  437. C  LIW -- Set it to the declared length of the IWORK array.
  438. C               You must have LIW .GE. 20+NEQ
  439. C
  440. C  RPAR, IPAR -- These are parameter arrays, of real and integer
  441. C         type, respectively.  You can use them for communication
  442. C         between your program that calls DDASSL and the
  443. C         RES subroutine (and the JAC subroutine).  They are not
  444. C         altered by DDASSL.  If you do not need RPAR or IPAR,
  445. C         ignore these parameters by treating them as dummy
  446. C         arguments.  If you do choose to use them, dimension
  447. C         them in your calling program and in RES (and in JAC)
  448. C         as arrays of appropriate length.
  449. C
  450. C  JAC -- If you have set INFO(5)=0, you can ignore this parameter
  451. C         by treating it as a dummy argument.  Otherwise, you must
  452. C         provide a subroutine of the form
  453. C             SUBROUTINE JAC(T,Y,YPRIME,PD,CJ,RPAR,IPAR)
  454. C         to define the matrix of partial derivatives
  455. C             PD=DG/DY+CJ*DG/DYPRIME
  456. C         CJ is a scalar which is input to JAC.
  457. C         For the given values of T,Y,YPRIME, the
  458. C         subroutine must evaluate the non-zero partial
  459. C         derivatives for each equation and each solution
  460. C         component, and store these values in the
  461. C         matrix PD.  The elements of PD are set to zero
  462. C         before each call to JAC so only non-zero elements
  463. C         need to be defined.
  464. C
  465. C         Subroutine JAC must not alter T,Y,(*),YPRIME(*), or CJ.
  466. C         You must declare the name JAC in an EXTERNAL statement in
  467. C         your program that calls DDASSL.  You must dimension Y,
  468. C         YPRIME and PD in JAC.
  469. C
  470. C         The way you must store the elements into the PD matrix
  471. C         depends on the structure of the matrix which you
  472. C         indicated by INFO(6).
  473. C               *** INFO(6)=0 -- Full (dense) matrix ***
  474. C                   Give PD a first dimension of NEQ.
  475. C                   When you evaluate the (non-zero) partial derivative
  476. C                   of equation I with respect to variable J, you must
  477. C                   store it in PD according to
  478. C                   PD(I,J) = "DG(I)/DY(J)+CJ*DG(I)/DYPRIME(J)"
  479. C               *** INFO(6)=1 -- Banded JACOBIAN with ML lower and MU
  480. C                   upper diagonal bands (refer to INFO(6) description
  481. C                   of ML and MU) ***
  482. C                   Give PD a first dimension of 2*ML+MU+1.
  483. C                   when you evaluate the (non-zero) partial derivative
  484. C                   of equation I with respect to variable J, you must
  485. C                   store it in PD according to
  486. C                   IROW = I - J + ML + MU + 1
  487. C                   PD(IROW,J) = "DG(I)/DY(J)+CJ*DG(I)/DYPRIME(J)"
  488. C
  489. C         RPAR and IPAR are real and integer parameter arrays
  490. C         which you can use for communication between your calling
  491. C         program and your JACOBIAN subroutine JAC. They are not
  492. C         altered by DDASSL. If you do not need RPAR or IPAR,
  493. C         ignore these parameters by treating them as dummy
  494. C         arguments. If you do choose to use them, dimension
  495. C         them in your calling program and in JAC as arrays of
  496. C         appropriate length.
  497. C
  498. C
  499. C  OPTIONALLY REPLACEABLE NORM ROUTINE:
  500. C
  501. C     DDASSL uses a weighted norm DDANRM to measure the size
  502. C     of vectors such as the estimated error in each step.
  503. C     A FUNCTION subprogram
  504. C       DOUBLE PRECISION FUNCTION DDANRM(NEQ,V,WT,RPAR,IPAR)
  505. C       DIMENSION V(NEQ),WT(NEQ)
  506. C     is used to define this norm. Here, V is the vector
  507. C     whose norm is to be computed, and WT is a vector of
  508. C     weights.  A DDANRM routine has been included with DDASSL
  509. C     which computes the weighted root-mean-square norm
  510. C     given by
  511. C       DDANRM=SQRT((1/NEQ)*SUM(V(I)/WT(I))**2)
  512. C     this norm is suitable for most problems. In some
  513. C     special cases, it may be more convenient and/or
  514. C     efficient to define your own norm by writing a function
  515. C     subprogram to be called instead of DDANRM. This should,
  516. C     however, be attempted only after careful thought and
  517. C     consideration.
  518. C
  519. C
  520. C  -------- OUTPUT -- AFTER ANY RETURN FROM DDASSL ---------------------
  521. C
  522. C  The principal aim of the code is to return a computed solution at
  523. C  TOUT, although it is also possible to obtain intermediate results
  524. C  along the way. To find out whether the code achieved its goal
  525. C  or if the integration process was interrupted before the task was
  526. C  completed, you must check the IDID parameter.
  527. C
  528. C
  529. C  T -- The solution was successfully advanced to the
  530. C               output value of T.
  531. C
  532. C  Y(*) -- Contains the computed solution approximation at T.
  533. C
  534. C  YPRIME(*) -- Contains the computed derivative
  535. C               approximation at T.
  536. C
  537. C  IDID -- Reports what the code did.
  538. C
  539. C                     *** Task completed ***
  540. C                Reported by positive values of IDID
  541. C
  542. C           IDID = 1 -- A step was successfully taken in the
  543. C                   intermediate-output mode. The code has not
  544. C                   yet reached TOUT.
  545. C
  546. C           IDID = 2 -- The integration to TSTOP was successfully
  547. C                   completed (T=TSTOP) by stepping exactly to TSTOP.
  548. C
  549. C           IDID = 3 -- The integration to TOUT was successfully
  550. C                   completed (T=TOUT) by stepping past TOUT.
  551. C                   Y(*) is obtained by interpolation.
  552. C                   YPRIME(*) is obtained by interpolation.
  553. C
  554. C                    *** Task interrupted ***
  555. C                Reported by negative values of IDID
  556. C
  557. C           IDID = -1 -- A large amount of work has been expended.
  558. C                   (About 500 steps)
  559. C
  560. C           IDID = -2 -- The error tolerances are too stringent.
  561. C
  562. C           IDID = -3 -- The local error test cannot be satisfied
  563. C                   because you specified a zero component in ATOL
  564. C                   and the corresponding computed solution
  565. C                   component is zero. Thus, a pure relative error
  566. C                   test is impossible for this component.
  567. C
  568. C           IDID = -6 -- DDASSL had repeated error test
  569. C                   failures on the last attempted step.
  570. C
  571. C           IDID = -7 -- The corrector could not converge.
  572. C
  573. C           IDID = -8 -- The matrix of partial derivatives
  574. C                   is singular.
  575. C
  576. C           IDID = -9 -- The corrector could not converge.
  577. C                   there were repeated error test failures
  578. C                   in this step.
  579. C
  580. C           IDID =-10 -- The corrector could not converge
  581. C                   because IRES was equal to minus one.
  582. C
  583. C           IDID =-11 -- IRES equal to -2 was encountered
  584. C                   and control is being returned to the
  585. C                   calling program.
  586. C
  587. C           IDID =-12 -- DDASSL failed to compute the initial
  588. C                   YPRIME.
  589. C
  590. C
  591. C
  592. C           IDID = -13,..,-32 -- Not applicable for this code
  593. C
  594. C                    *** Task terminated ***
  595. C                Reported by the value of IDID=-33
  596. C
  597. C           IDID = -33 -- The code has encountered trouble from which
  598. C                   it cannot recover. A message is printed
  599. C                   explaining the trouble and control is returned
  600. C                   to the calling program. For example, this occurs
  601. C                   when invalid input is detected.
  602. C
  603. C  RTOL, ATOL -- These quantities remain unchanged except when
  604. C               IDID = -2. In this case, the error tolerances have been
  605. C               increased by the code to values which are estimated to
  606. C               be appropriate for continuing the integration. However,
  607. C               the reported solution at T was obtained using the input
  608. C               values of RTOL and ATOL.
  609. C
  610. C  RWORK, IWORK -- Contain information which is usually of no
  611. C               interest to the user but necessary for subsequent calls.
  612. C               However, you may find use for
  613. C
  614. C               RWORK(3)--Which contains the step size H to be
  615. C                       attempted on the next step.
  616. C
  617. C               RWORK(4)--Which contains the current value of the
  618. C                       independent variable, i.e., the farthest point
  619. C                       integration has reached. This will be different
  620. C                       from T only when interpolation has been
  621. C                       performed (IDID=3).
  622. C
  623. C               RWORK(7)--Which contains the stepsize used
  624. C                       on the last successful step.
  625. C
  626. C               IWORK(7)--Which contains the order of the method to
  627. C                       be attempted on the next step.
  628. C
  629. C               IWORK(8)--Which contains the order of the method used
  630. C                       on the last step.
  631. C
  632. C               IWORK(11)--Which contains the number of steps taken so
  633. C                        far.
  634. C
  635. C               IWORK(12)--Which contains the number of calls to RES
  636. C                        so far.
  637. C
  638. C               IWORK(13)--Which contains the number of evaluations of
  639. C                        the matrix of partial derivatives needed so
  640. C                        far.
  641. C
  642. C               IWORK(14)--Which contains the total number
  643. C                        of error test failures so far.
  644. C
  645. C               IWORK(15)--Which contains the total number
  646. C                        of convergence test failures so far.
  647. C                        (includes singular iteration matrix
  648. C                        failures.)
  649. C
  650. C
  651. C  -------- INPUT -- WHAT TO DO TO CONTINUE THE INTEGRATION ------------
  652. C                    (CALLS AFTER THE FIRST)
  653. C
  654. C  This code is organized so that subsequent calls to continue the
  655. C  integration involve little (if any) additional effort on your
  656. C  part. You must monitor the IDID parameter in order to determine
  657. C  what to do next.
  658. C
  659. C  Recalling that the principal task of the code is to integrate
  660. C  from T to TOUT (the interval mode), usually all you will need
  661. C  to do is specify a new TOUT upon reaching the current TOUT.
  662. C
  663. C  Do not alter any quantity not specifically permitted below,
  664. C  in particular do not alter NEQ,T,Y(*),YPRIME(*),RWORK(*),IWORK(*)
  665. C  or the differential equation in subroutine RES. Any such
  666. C  alteration constitutes a new problem and must be treated as such,
  667. C  i.e., you must start afresh.
  668. C
  669. C  You cannot change from vector to scalar error control or vice
  670. C  versa (INFO(2)), but you can change the size of the entries of
  671. C  RTOL, ATOL. Increasing a tolerance makes the equation easier
  672. C  to integrate. Decreasing a tolerance will make the equation
  673. C  harder to integrate and should generally be avoided.
  674. C
  675. C  You can switch from the intermediate-output mode to the
  676. C  interval mode (INFO(3)) or vice versa at any time.
  677. C
  678. C  If it has been necessary to prevent the integration from going
  679. C  past a point TSTOP (INFO(4), RWORK(1)), keep in mind that the
  680. C  code will not integrate to any TOUT beyond the currently
  681. C  specified TSTOP. Once TSTOP has been reached you must change
  682. C  the value of TSTOP or set INFO(4)=0. You may change INFO(4)
  683. C  or TSTOP at any time but you must supply the value of TSTOP in
  684. C  RWORK(1) whenever you set INFO(4)=1.
  685. C
  686. C  Do not change INFO(5), INFO(6), IWORK(1), or IWORK(2)
  687. C  unless you are going to restart the code.
  688. C
  689. C                 *** Following a completed task ***
  690. C  If
  691. C     IDID = 1, call the code again to continue the integration
  692. C                  another step in the direction of TOUT.
  693. C
  694. C     IDID = 2 or 3, define a new TOUT and call the code again.
  695. C                  TOUT must be different from T. You cannot change
  696. C                  the direction of integration without restarting.
  697. C
  698. C                 *** Following an interrupted task ***
  699. C               To show the code that you realize the task was
  700. C               interrupted and that you want to continue, you
  701. C               must take appropriate action and set INFO(1) = 1
  702. C  If
  703. C    IDID = -1, The code has taken about 500 steps.
  704. C                  If you want to continue, set INFO(1) = 1 and
  705. C                  call the code again. An additional 500 steps
  706. C                  will be allowed.
  707. C
  708. C    IDID = -2, The error tolerances RTOL, ATOL have been
  709. C                  increased to values the code estimates appropriate
  710. C                  for continuing. You may want to change them
  711. C                  yourself. If you are sure you want to continue
  712. C                  with relaxed error tolerances, set INFO(1)=1 and
  713. C                  call the code again.
  714. C
  715. C    IDID = -3, A solution component is zero and you set the
  716. C                  corresponding component of ATOL to zero. If you
  717. C                  are sure you want to continue, you must first
  718. C                  alter the error criterion to use positive values
  719. C                  for those components of ATOL corresponding to zero
  720. C                  solution components, then set INFO(1)=1 and call
  721. C                  the code again.
  722. C
  723. C    IDID = -4,-5  --- Cannot occur with this code.
  724. C
  725. C    IDID = -6, Repeated error test failures occurred on the
  726. C                  last attempted step in DDASSL. A singularity in the
  727. C                  solution may be present. If you are absolutely
  728. C                  certain you want to continue, you should restart
  729. C                  the integration. (Provide initial values of Y and
  730. C                  YPRIME which are consistent)
  731. C
  732. C    IDID = -7, Repeated convergence test failures occurred
  733. C                  on the last attempted step in DDASSL. An inaccurate
  734. C                  or ill-conditioned JACOBIAN may be the problem. If
  735. C                  you are absolutely certain you want to continue, you
  736. C                  should restart the integration.
  737. C
  738. C    IDID = -8, The matrix of partial derivatives is singular.
  739. C                  Some of your equations may be redundant.
  740. C                  DDASSL cannot solve the problem as stated.
  741. C                  It is possible that the redundant equations
  742. C                  could be removed, and then DDASSL could
  743. C                  solve the problem. It is also possible
  744. C                  that a solution to your problem either
  745. C                  does not exist or is not unique.
  746. C
  747. C    IDID = -9, DDASSL had multiple convergence test
  748. C                  failures, preceeded by multiple error
  749. C                  test failures, on the last attempted step.
  750. C                  It is possible that your problem
  751. C                  is ill-posed, and cannot be solved
  752. C                  using this code. Or, there may be a
  753. C                  discontinuity or a singularity in the
  754. C                  solution. If you are absolutely certain
  755. C                  you want to continue, you should restart
  756. C                  the integration.
  757. C
  758. C    IDID =-10, DDASSL had multiple convergence test failures
  759. C                  because IRES was equal to minus one.
  760. C                  If you are absolutely certain you want
  761. C                  to continue, you should restart the
  762. C                  integration.
  763. C
  764. C    IDID =-11, IRES=-2 was encountered, and control is being
  765. C                  returned to the calling program.
  766. C
  767. C    IDID =-12, DDASSL failed to compute the initial YPRIME.
  768. C                  This could happen because the initial
  769. C                  approximation to YPRIME was not very good, or
  770. C                  if a YPRIME consistent with the initial Y
  771. C                  does not exist. The problem could also be caused
  772. C                  by an inaccurate or singular iteration matrix.
  773. C
  774. C    IDID = -13,..,-32  --- Cannot occur with this code.
  775. C
  776. C
  777. C                 *** Following a terminated task ***
  778. C
  779. C  If IDID= -33, you cannot continue the solution of this problem.
  780. C                  An attempt to do so will result in your
  781. C                  run being terminated.
  782. C
  783. C
  784. C  -------- ERROR MESSAGES ---------------------------------------------
  785. C
  786. C      The SLATEC error print routine XERMSG is called in the event of
  787. C   unsuccessful completion of a task.  Most of these are treated as
  788. C   "recoverable errors", which means that (unless the user has directed
  789. C   otherwise) control will be returned to the calling program for
  790. C   possible action after the message has been printed.
  791. C
  792. C   In the event of a negative value of IDID other than -33, an appro-
  793. C   priate message is printed and the "error number" printed by XERMSG
  794. C   is the value of IDID.  There are quite a number of illegal input
  795. C   errors that can lead to a returned value IDID=-33.  The conditions
  796. C   and their printed "error numbers" are as follows:
  797. C
  798. C   Error number       Condition
  799. C
  800. C        1       Some element of INFO vector is not zero or one.
  801. C        2       NEQ .le. 0
  802. C        3       MAXORD not in range.
  803. C        4       LRW is less than the required length for RWORK.
  804. C        5       LIW is less than the required length for IWORK.
  805. C        6       Some element of RTOL is .lt. 0
  806. C        7       Some element of ATOL is .lt. 0
  807. C        8       All elements of RTOL and ATOL are zero.
  808. C        9       INFO(4)=1 and TSTOP is behind TOUT.
  809. C       10       HMAX .lt. 0.0
  810. C       11       TOUT is behind T.
  811. C       12       INFO(8)=1 and H0=0.0
  812. C       13       Some element of WT is .le. 0.0
  813. C       14       TOUT is too close to T to start integration.
  814. C       15       INFO(4)=1 and TSTOP is behind T.
  815. C       16       --( Not used in this version )--
  816. C       17       ML illegal.  Either .lt. 0 or .gt. NEQ
  817. C       18       MU illegal.  Either .lt. 0 or .gt. NEQ
  818. C       19       TOUT = T.
  819. C
  820. C   If DDASSL is called again without any action taken to remove the
  821. C   cause of an unsuccessful return, XERMSG will be called with a fatal
  822. C   error flag, which will cause unconditional termination of the
  823. C   program.  There are two such fatal errors:
  824. C
  825. C   Error number -998:  The last step was terminated with a negative
  826. C       value of IDID other than -33, and no appropriate action was
  827. C       taken.
  828. C
  829. C   Error number -999:  The previous call was terminated because of
  830. C       illegal input (IDID=-33) and there is illegal input in the
  831. C       present call, as well.  (Suspect infinite loop.)
  832. C
  833. C  ---------------------------------------------------------------------
  834. C
  835. C***REFERENCES  A DESCRIPTION OF DASSL: A DIFFERENTIAL/ALGEBRAIC
  836. C                 SYSTEM SOLVER, L. R. PETZOLD, SAND82-8637,
  837. C                 SANDIA NATIONAL LABORATORIES, SEPTEMBER 1982.
  838. C***ROUTINES CALLED  D1MACH, DDAINI, DDANRM, DDASTP, DDATRP, DDAWTS,
  839. C                    XERMSG
  840. C***REVISION HISTORY  (YYMMDD)
  841. C   830315  DATE WRITTEN
  842. C   880387  Code changes made.  All common statements have been
  843. C           replaced by a DATA statement, which defines pointers into
  844. C           RWORK, and PARAMETER statements which define pointers
  845. C           into IWORK.  As well the documentation has gone through
  846. C           grammatical changes.
  847. C   881005  The prologue has been changed to mixed case.
  848. C           The subordinate routines had revision dates changed to
  849. C           this date, although the documentation for these routines
  850. C           is all upper case.  No code changes.
  851. C   890511  Code changes made.  The DATA statement in the declaration
  852. C           section of DDASSL was replaced with a PARAMETER
  853. C           statement.  Also the statement S = 100.D0 was removed
  854. C           from the top of the Newton iteration in DDASTP.
  855. C           The subordinate routines had revision dates changed to
  856. C           this date.
  857. C   890517  The revision date syntax was replaced with the revision
  858. C           history syntax.  Also the "DECK" comment was added to
  859. C           the top of all subroutines.  These changes are consistent
  860. C           with new SLATEC guidelines.
  861. C           The subordinate routines had revision dates changed to
  862. C           this date.  No code changes.
  863. C   891013  Code changes made.
  864. C           Removed all occurrances of FLOAT or DBLE.  All operations
  865. C           are now performed with "mixed-mode" arithmetic.
  866. C           Also, specific function names were replaced with generic
  867. C           function names to be consistent with new SLATEC guidelines.
  868. C           In particular:
  869. C              Replaced DSQRT with SQRT everywhere.
  870. C              Replaced DABS with ABS everywhere.
  871. C              Replaced DMIN1 with MIN everywhere.
  872. C              Replaced MIN0 with MIN everywhere.
  873. C              Replaced DMAX1 with MAX everywhere.
  874. C              Replaced MAX0 with MAX everywhere.
  875. C              Replaced DSIGN with SIGN everywhere.
  876. C           Also replaced REVISION DATE with REVISION HISTORY in all
  877. C           subordinate routines.
  878. C  901004  Miscellaneous changes to prologue to complete conversion
  879. C          to SLATEC 4.0 format.  No code changes.  (F.N.Fritsch)
  880. C  901009  Corrected GAMS classification code and converted subsidiary
  881. C          routines to 4.0 format.  No code changes.  (F.N.Fritsch)
  882. C  901010  Converted XERRWV calls to XERMSG calls.  (R.Clemens,AFWL)
  883. C  901019  Code changes made.
  884. C          Merged SLATEC 4.0 changes with previous changes made
  885. C          by C. Ulrich.  Below is a history of the changes made by
  886. C          C. Ulrich. (Changes in subsidiary routines are implied
  887. C          by this history)
  888. C          891228  Bug was found and repaired inside the DDASSL
  889. C                  and DDAINI routines.  DDAINI was incorrectly
  890. C                  returning the initial T with Y and YPRIME
  891. C                  computed at T+H.  The routine now returns T+H
  892. C                  rather than the initial T.
  893. C                  Cosmetic changes made to DDASTP.
  894. C          900904  Three modifications were made to fix a bug (inside
  895. C                  DDASSL) re interpolation for continuation calls and
  896. C                  cases where TN is very close to TSTOP:
  897. C
  898. C                  1) In testing for whether H is too large, just
  899. C                     compare H to (TSTOP - TN), rather than
  900. C                     (TSTOP - TN) * (1-4*UROUND), and set H to
  901. C                     TSTOP - TN.  This will force DDASTP to step
  902. C                     exactly to TSTOP under certain situations
  903. C                     (i.e. when H returned from DDASTP would otherwise
  904. C                     take TN beyond TSTOP).
  905. C
  906. C                  2) Inside the DDASTP loop, interpolate exactly to
  907. C                     TSTOP if TN is very close to TSTOP (rather than
  908. C                     interpolating to within roundoff of TSTOP).
  909. C
  910. C                  3) Modified IDID description for IDID = 2 to say that
  911. C                     the solution is returned by stepping exactly to
  912. C                     TSTOP, rather than TOUT.  (In some cases the
  913. C                     solution is actually obtained by extrapolating
  914. C                     over a distance near unit roundoff to TSTOP,
  915. C                     but this small distance is deemed acceptable in
  916. C                     these circumstances.)
  917. C   901026  Added explicit declarations for all variables and minor
  918. C           cosmetic changes to prologue, removed unreferenced labels,
  919. C           and improved XERMSG calls.  (FNF)
  920. C   901030  Added ERROR MESSAGES section and reworked other sections to
  921. C           be of more uniform format.  (FNF)
  922. C   910624  Fixed minor bug related to HMAX (five lines ending in
  923. C           statement 526 in DDASSL).   (LRP)
  924. C
  925. C***END PROLOGUE  DDASSL
  926. C
  927. C**End
  928. C
  929. C     Declare arguments.
  930. C
  931.       INTEGER  NEQ, INFO(15), IDID, LRW, IWORK(*), LIW, IPAR(*)
  932.       DOUBLE PRECISION
  933.      *   T, Y(*), YPRIME(*), TOUT, RTOL(*), ATOL(*), RWORK(*),
  934.      *   RPAR(*)
  935.       EXTERNAL  RES, JAC
  936. C
  937. C     Declare externals.
  938. C
  939.       EXTERNAL  D1MACH, DDAINI, DDANRM, DDASTP, DDATRP, DDAWTS, XERMSG
  940.       DOUBLE PRECISION  D1MACH, DDANRM
  941. C
  942. C     Declare local variables.
  943. C
  944.       INTEGER  I, ITEMP, LALPHA, LBETA, LCJ, LCJOLD, LCTF, LDELTA,
  945.      *   LENIW, LENPD, LENRW, LE, LETF, LGAMMA, LH, LHMAX, LHOLD, LIPVT,
  946.      *   LJCALC, LK, LKOLD, LIWM, LML, LMTYPE, LMU, LMXORD, LNJE, LNPD,
  947.      *   LNRE, LNS, LNST, LNSTL, LPD, LPHASE, LPHI, LPSI, LROUND, LS,
  948.      *   LSIGMA, LTN, LTSTOP, LWM, LWT, MBAND, MSAVE, MXORD, NPD, NTEMP,
  949.      *   NZFLG
  950.       DOUBLE PRECISION
  951.      *   ATOLI, H, HMAX, HMIN, HO, R, RH, RTOLI, TDIST, TN, TNEXT,
  952.      *   TSTOP, UROUND, YPNORM
  953.       LOGICAL  DONE
  954. C       Auxiliary variables for conversion of values to be included in
  955. C       error messages.
  956.       CHARACTER*8  XERN1, XERN2
  957.       CHARACTER*16 XERN3, XERN4
  958. C
  959. C     SET POINTERS INTO IWORK
  960.       PARAMETER (LML=1, LMU=2, LMXORD=3, LMTYPE=4, LNST=11,
  961.      *  LNRE=12, LNJE=13, LETF=14, LCTF=15, LNPD=16,
  962.      *  LIPVT=21, LJCALC=5, LPHASE=6, LK=7, LKOLD=8,
  963.      *  LNS=9, LNSTL=10, LIWM=1)
  964. C
  965. C     SET RELATIVE OFFSET INTO RWORK
  966.       PARAMETER (NPD=1)
  967. C
  968. C     SET POINTERS INTO RWORK
  969.       PARAMETER (LTSTOP=1, LHMAX=2, LH=3, LTN=4,
  970.      *  LCJ=5, LCJOLD=6, LHOLD=7, LS=8, LROUND=9,
  971.      *  LALPHA=11, LBETA=17, LGAMMA=23,
  972.      *  LPSI=29, LSIGMA=35, LDELTA=41)
  973. C
  974. C***FIRST EXECUTABLE STATEMENT  DDASSL
  975.       IF(INFO(1).NE.0)GO TO 100
  976. C
  977. C-----------------------------------------------------------------------
  978. C     THIS BLOCK IS EXECUTED FOR THE INITIAL CALL ONLY.
  979. C     IT CONTAINS CHECKING OF INPUTS AND INITIALIZATIONS.
  980. C-----------------------------------------------------------------------
  981. C
  982. C     FIRST CHECK INFO ARRAY TO MAKE SURE ALL ELEMENTS OF INFO
  983. C     ARE EITHER ZERO OR ONE.
  984.       DO 10 I=2,11
  985.          IF(INFO(I).NE.0.AND.INFO(I).NE.1)GO TO 701
  986. 10       CONTINUE
  987. C
  988.       IF(NEQ.LE.0)GO TO 702
  989. C
  990. C     CHECK AND COMPUTE MAXIMUM ORDER
  991.       MXORD=5
  992.       IF(INFO(9).EQ.0)GO TO 20
  993.          MXORD=IWORK(LMXORD)
  994.          IF(MXORD.LT.1.OR.MXORD.GT.5)GO TO 703
  995. 20       IWORK(LMXORD)=MXORD
  996. C
  997. C     COMPUTE MTYPE,LENPD,LENRW.CHECK ML AND MU.
  998.       IF(INFO(6).NE.0)GO TO 40
  999.          LENPD=NEQ**2
  1000.          LENRW=40+(IWORK(LMXORD)+4)*NEQ+LENPD
  1001.          IF(INFO(5).NE.0)GO TO 30
  1002.             IWORK(LMTYPE)=2
  1003.             GO TO 60
  1004. 30          IWORK(LMTYPE)=1
  1005.             GO TO 60
  1006. 40    IF(IWORK(LML).LT.0.OR.IWORK(LML).GE.NEQ)GO TO 717
  1007.       IF(IWORK(LMU).LT.0.OR.IWORK(LMU).GE.NEQ)GO TO 718
  1008.       LENPD=(2*IWORK(LML)+IWORK(LMU)+1)*NEQ
  1009.       IF(INFO(5).NE.0)GO TO 50
  1010.          IWORK(LMTYPE)=5
  1011.          MBAND=IWORK(LML)+IWORK(LMU)+1
  1012.          MSAVE=(NEQ/MBAND)+1
  1013.          LENRW=40+(IWORK(LMXORD)+4)*NEQ+LENPD+2*MSAVE
  1014.          GO TO 60
  1015. 50       IWORK(LMTYPE)=4
  1016.          LENRW=40+(IWORK(LMXORD)+4)*NEQ+LENPD
  1017. C
  1018. C     CHECK LENGTHS OF RWORK AND IWORK
  1019. 60    LENIW=20+NEQ
  1020.       IWORK(LNPD)=LENPD
  1021.       IF(LRW.LT.LENRW)GO TO 704
  1022.       IF(LIW.LT.LENIW)GO TO 705
  1023. C
  1024. C     CHECK TO SEE THAT TOUT IS DIFFERENT FROM T
  1025.       IF(TOUT .EQ. T)GO TO 719
  1026. C
  1027. C     CHECK HMAX
  1028.       IF(INFO(7).EQ.0)GO TO 70
  1029.          HMAX=RWORK(LHMAX)
  1030.          IF(HMAX.LE.0.0D0)GO TO 710
  1031. 70    CONTINUE
  1032. C
  1033. C     INITIALIZE COUNTERS
  1034.       IWORK(LNST)=0
  1035.       IWORK(LNRE)=0
  1036.       IWORK(LNJE)=0
  1037. C
  1038.       IWORK(LNSTL)=0
  1039.       IDID=1
  1040.       GO TO 200
  1041. C
  1042. C-----------------------------------------------------------------------
  1043. C     THIS BLOCK IS FOR CONTINUATION CALLS
  1044. C     ONLY. HERE WE CHECK INFO(1),AND IF THE
  1045. C     LAST STEP WAS INTERRUPTED WE CHECK WHETHER
  1046. C     APPROPRIATE ACTION WAS TAKEN.
  1047. C-----------------------------------------------------------------------
  1048. C
  1049. 100   CONTINUE
  1050.       IF(INFO(1).EQ.1)GO TO 110
  1051.       IF(INFO(1).NE.-1)GO TO 701
  1052. C
  1053. C     IF WE ARE HERE, THE LAST STEP WAS INTERRUPTED
  1054. C     BY AN ERROR CONDITION FROM DDASTP,AND
  1055. C     APPROPRIATE ACTION WAS NOT TAKEN. THIS
  1056. C     IS A FATAL ERROR.
  1057.       WRITE (XERN1, '(I8)') IDID
  1058.       CALL XERMSG ('SLATEC', 'DDASSL',
  1059.      *   'THE LAST STEP TERMINATED WITH A NEGATIVE VALUE OF IDID = ' //
  1060.      *   XERN1 // ' AND NO APPROPRIATE ACTION WAS TAKEN.  ' //
  1061.      *   'RUN TERMINATED', -998, 2)
  1062.       RETURN
  1063. 110   CONTINUE
  1064.       IWORK(LNSTL)=IWORK(LNST)
  1065. C
  1066. C-----------------------------------------------------------------------
  1067. C     THIS BLOCK IS EXECUTED ON ALL CALLS.
  1068. C     THE ERROR TOLERANCE PARAMETERS ARE
  1069. C     CHECKED, AND THE WORK ARRAY POINTERS
  1070. C     ARE SET.
  1071. C-----------------------------------------------------------------------
  1072. C
  1073. 200   CONTINUE
  1074. C     CHECK RTOL,ATOL
  1075.       NZFLG=0
  1076.       RTOLI=RTOL(1)
  1077.       ATOLI=ATOL(1)
  1078.       DO 210 I=1,NEQ
  1079.          IF(INFO(2).EQ.1)RTOLI=RTOL(I)
  1080.          IF(INFO(2).EQ.1)ATOLI=ATOL(I)
  1081.          IF(RTOLI.GT.0.0D0.OR.ATOLI.GT.0.0D0)NZFLG=1
  1082.          IF(RTOLI.LT.0.0D0)GO TO 706
  1083.          IF(ATOLI.LT.0.0D0)GO TO 707
  1084. 210      CONTINUE
  1085.       IF(NZFLG.EQ.0)GO TO 708
  1086. C
  1087. C     SET UP RWORK STORAGE.IWORK STORAGE IS FIXED
  1088. C     IN DATA STATEMENT.
  1089.       LE=LDELTA+NEQ
  1090.       LWT=LE+NEQ
  1091.       LPHI=LWT+NEQ
  1092.       LPD=LPHI+(IWORK(LMXORD)+1)*NEQ
  1093.       LWM=LPD
  1094.       NTEMP=NPD+IWORK(LNPD)
  1095.       IF(INFO(1).EQ.1)GO TO 400
  1096. C
  1097. C-----------------------------------------------------------------------
  1098. C     THIS BLOCK IS EXECUTED ON THE INITIAL CALL
  1099. C     ONLY. SET THE INITIAL STEP SIZE, AND
  1100. C     THE ERROR WEIGHT VECTOR, AND PHI.
  1101. C     COMPUTE INITIAL YPRIME, IF NECESSARY.
  1102. C-----------------------------------------------------------------------
  1103. C
  1104.       TN=T
  1105.       IDID=1
  1106. C
  1107. C     SET ERROR WEIGHT VECTOR WT
  1108.       CALL DDAWTS(NEQ,INFO(2),RTOL,ATOL,Y,RWORK(LWT),RPAR,IPAR)
  1109.       DO 305 I = 1,NEQ
  1110.          IF(RWORK(LWT+I-1).LE.0.0D0) GO TO 713
  1111. 305      CONTINUE
  1112. C
  1113. C     COMPUTE UNIT ROUNDOFF AND HMIN
  1114.       UROUND = D1MACH(4)
  1115.       RWORK(LROUND) = UROUND
  1116.       HMIN = 4.0D0*UROUND*MAX(ABS(T),ABS(TOUT))
  1117. C
  1118. C     CHECK INITIAL INTERVAL TO SEE THAT IT IS LONG ENOUGH
  1119.       TDIST = ABS(TOUT - T)
  1120.       IF(TDIST .LT. HMIN) GO TO 714
  1121. C
  1122. C     CHECK HO, IF THIS WAS INPUT
  1123.       IF (INFO(8) .EQ. 0) GO TO 310
  1124.          HO = RWORK(LH)
  1125.          IF ((TOUT - T)*HO .LT. 0.0D0) GO TO 711
  1126.          IF (HO .EQ. 0.0D0) GO TO 712
  1127.          GO TO 320
  1128. 310    CONTINUE
  1129. C
  1130. C     COMPUTE INITIAL STEPSIZE, TO BE USED BY EITHER
  1131. C     DDASTP OR DDAINI, DEPENDING ON INFO(11)
  1132.       HO = 0.001D0*TDIST
  1133.       YPNORM = DDANRM(NEQ,YPRIME,RWORK(LWT),RPAR,IPAR)
  1134.       IF (YPNORM .GT. 0.5D0/HO) HO = 0.5D0/YPNORM
  1135.       HO = SIGN(HO,TOUT-T)
  1136. C     ADJUST HO IF NECESSARY TO MEET HMAX BOUND
  1137. 320   IF (INFO(7) .EQ. 0) GO TO 330
  1138.          RH = ABS(HO)/RWORK(LHMAX)
  1139.          IF (RH .GT. 1.0D0) HO = HO/RH
  1140. C     COMPUTE TSTOP, IF APPLICABLE
  1141. 330   IF (INFO(4) .EQ. 0) GO TO 340
  1142.          TSTOP = RWORK(LTSTOP)
  1143.          IF ((TSTOP - T)*HO .LT. 0.0D0) GO TO 715
  1144.          IF ((T + HO - TSTOP)*HO .GT. 0.0D0) HO = TSTOP - T
  1145.          IF ((TSTOP - TOUT)*HO .LT. 0.0D0) GO TO 709
  1146. C
  1147. C     COMPUTE INITIAL DERIVATIVE, UPDATING TN AND Y, IF APPLICABLE
  1148. 340   IF (INFO(11) .EQ. 0) GO TO 350
  1149.       CALL DDAINI(TN,Y,YPRIME,NEQ,
  1150.      *  RES,JAC,HO,RWORK(LWT),IDID,RPAR,IPAR,
  1151.      *  RWORK(LPHI),RWORK(LDELTA),RWORK(LE),
  1152.      *  RWORK(LWM),IWORK(LIWM),HMIN,RWORK(LROUND),
  1153.      *  INFO(10),NTEMP)
  1154.       IF (IDID .LT. 0) GO TO 390
  1155. C
  1156. C     LOAD H WITH HO.  STORE H IN RWORK(LH)
  1157. 350   H = HO
  1158.       RWORK(LH) = H
  1159. C
  1160. C     LOAD Y AND H*YPRIME INTO PHI(*,1) AND PHI(*,2)
  1161.       ITEMP = LPHI + NEQ
  1162.       DO 370 I = 1,NEQ
  1163.          RWORK(LPHI + I - 1) = Y(I)
  1164. 370      RWORK(ITEMP + I - 1) = H*YPRIME(I)
  1165. C
  1166. 390   GO TO 500
  1167. C
  1168. C-------------------------------------------------------
  1169. C     THIS BLOCK IS FOR CONTINUATION CALLS ONLY. ITS
  1170. C     PURPOSE IS TO CHECK STOP CONDITIONS BEFORE
  1171. C     TAKING A STEP.
  1172. C     ADJUST H IF NECESSARY TO MEET HMAX BOUND
  1173. C-------------------------------------------------------
  1174. C
  1175. 400   CONTINUE
  1176.       UROUND=RWORK(LROUND)
  1177.       DONE = .FALSE.
  1178.       TN=RWORK(LTN)
  1179.       H=RWORK(LH)
  1180.       IF(INFO(7) .EQ. 0) GO TO 410
  1181.          RH = ABS(H)/RWORK(LHMAX)
  1182.          IF(RH .GT. 1.0D0) H = H/RH
  1183. 410   CONTINUE
  1184.       IF(T .EQ. TOUT) GO TO 719
  1185.       IF((T - TOUT)*H .GT. 0.0D0) GO TO 711
  1186.       IF(INFO(4) .EQ. 1) GO TO 430
  1187.       IF(INFO(3) .EQ. 1) GO TO 420
  1188.       IF((TN-TOUT)*H.LT.0.0D0)GO TO 490
  1189.       CALL DDATRP(TN,TOUT,Y,YPRIME,NEQ,IWORK(LKOLD),
  1190.      *  RWORK(LPHI),RWORK(LPSI))
  1191.       T=TOUT
  1192.       IDID = 3
  1193.       DONE = .TRUE.
  1194.       GO TO 490
  1195. 420   IF((TN-T)*H .LE. 0.0D0) GO TO 490
  1196.       IF((TN - TOUT)*H .GT. 0.0D0) GO TO 425
  1197.       CALL DDATRP(TN,TN,Y,YPRIME,NEQ,IWORK(LKOLD),
  1198.      *  RWORK(LPHI),RWORK(LPSI))
  1199.       T = TN
  1200.       IDID = 1
  1201.       DONE = .TRUE.
  1202.       GO TO 490
  1203. 425   CONTINUE
  1204.       CALL DDATRP(TN,TOUT,Y,YPRIME,NEQ,IWORK(LKOLD),
  1205.      *  RWORK(LPHI),RWORK(LPSI))
  1206.       T = TOUT
  1207.       IDID = 3
  1208.       DONE = .TRUE.
  1209.       GO TO 490
  1210. 430   IF(INFO(3) .EQ. 1) GO TO 440
  1211.       TSTOP=RWORK(LTSTOP)
  1212.       IF((TN-TSTOP)*H.GT.0.0D0) GO TO 715
  1213.       IF((TSTOP-TOUT)*H.LT.0.0D0)GO TO 709
  1214.       IF((TN-TOUT)*H.LT.0.0D0)GO TO 450
  1215.       CALL DDATRP(TN,TOUT,Y,YPRIME,NEQ,IWORK(LKOLD),
  1216.      *   RWORK(LPHI),RWORK(LPSI))
  1217.       T=TOUT
  1218.       IDID = 3
  1219.       DONE = .TRUE.
  1220.       GO TO 490
  1221. 440   TSTOP = RWORK(LTSTOP)
  1222.       IF((TN-TSTOP)*H .GT. 0.0D0) GO TO 715
  1223.       IF((TSTOP-TOUT)*H .LT. 0.0D0) GO TO 709
  1224.       IF((TN-T)*H .LE. 0.0D0) GO TO 450
  1225.       IF((TN - TOUT)*H .GT. 0.0D0) GO TO 445
  1226.       CALL DDATRP(TN,TN,Y,YPRIME,NEQ,IWORK(LKOLD),
  1227.      *  RWORK(LPHI),RWORK(LPSI))
  1228.       T = TN
  1229.       IDID = 1
  1230.       DONE = .TRUE.
  1231.       GO TO 490
  1232. 445   CONTINUE
  1233.       CALL DDATRP(TN,TOUT,Y,YPRIME,NEQ,IWORK(LKOLD),
  1234.      *  RWORK(LPHI),RWORK(LPSI))
  1235.       T = TOUT
  1236.       IDID = 3
  1237.       DONE = .TRUE.
  1238.       GO TO 490
  1239. 450   CONTINUE
  1240. C     CHECK WHETHER WE ARE WITHIN ROUNDOFF OF TSTOP
  1241.       IF(ABS(TN-TSTOP).GT.100.0D0*UROUND*
  1242.      *   (ABS(TN)+ABS(H)))GO TO 460
  1243.       CALL DDATRP(TN,TSTOP,Y,YPRIME,NEQ,IWORK(LKOLD),
  1244.      *  RWORK(LPHI),RWORK(LPSI))
  1245.       IDID=2
  1246.       T=TSTOP
  1247.       DONE = .TRUE.
  1248.       GO TO 490
  1249. 460   TNEXT=TN+H
  1250.       IF((TNEXT-TSTOP)*H.LE.0.0D0)GO TO 490
  1251.       H=TSTOP-TN
  1252.       RWORK(LH)=H
  1253. C
  1254. 490   IF (DONE) GO TO 580
  1255. C
  1256. C-------------------------------------------------------
  1257. C     THE NEXT BLOCK CONTAINS THE CALL TO THE
  1258. C     ONE-STEP INTEGRATOR DDASTP.
  1259. C     THIS IS A LOOPING POINT FOR THE INTEGRATION STEPS.
  1260. C     CHECK FOR TOO MANY STEPS.
  1261. C     UPDATE WT.
  1262. C     CHECK FOR TOO MUCH ACCURACY REQUESTED.
  1263. C     COMPUTE MINIMUM STEPSIZE.
  1264. C-------------------------------------------------------
  1265. C
  1266. 500   CONTINUE
  1267. C     CHECK FOR FAILURE TO COMPUTE INITIAL YPRIME
  1268.       IF (IDID .EQ. -12) GO TO 527
  1269. C
  1270. C     CHECK FOR TOO MANY STEPS
  1271.       IF((IWORK(LNST)-IWORK(LNSTL)).LT.500)
  1272.      *   GO TO 510
  1273.            IDID=-1
  1274.            GO TO 527
  1275. C
  1276. C     UPDATE WT
  1277. 510   CALL DDAWTS(NEQ,INFO(2),RTOL,ATOL,RWORK(LPHI),
  1278.      *  RWORK(LWT),RPAR,IPAR)
  1279.       DO 520 I=1,NEQ
  1280.          IF(RWORK(I+LWT-1).GT.0.0D0)GO TO 520
  1281.            IDID=-3
  1282.            GO TO 527
  1283. 520   CONTINUE
  1284. C
  1285. C     TEST FOR TOO MUCH ACCURACY REQUESTED.
  1286.       R=DDANRM(NEQ,RWORK(LPHI),RWORK(LWT),RPAR,IPAR)*
  1287.      *   100.0D0*UROUND
  1288.       IF(R.LE.1.0D0)GO TO 525
  1289. C     MULTIPLY RTOL AND ATOL BY R AND RETURN
  1290.       IF(INFO(2).EQ.1)GO TO 523
  1291.            RTOL(1)=R*RTOL(1)
  1292.            ATOL(1)=R*ATOL(1)
  1293.            IDID=-2
  1294.            GO TO 527
  1295. 523   DO 524 I=1,NEQ
  1296.            RTOL(I)=R*RTOL(I)
  1297. 524        ATOL(I)=R*ATOL(I)
  1298.       IDID=-2
  1299.       GO TO 527
  1300. 525   CONTINUE
  1301. C
  1302. C     COMPUTE MINIMUM STEPSIZE
  1303.       HMIN=4.0D0*UROUND*MAX(ABS(TN),ABS(TOUT))
  1304. C
  1305. C     TEST H VS. HMAX
  1306.       IF (INFO(7) .EQ. 0) GO TO 526
  1307.          RH = ABS(H)/RWORK(LHMAX)
  1308.          IF (RH .GT. 1.0D0) H = H/RH
  1309. 526   CONTINUE           
  1310. C
  1311.       CALL DDASTP(TN,Y,YPRIME,NEQ,
  1312.      *   RES,JAC,H,RWORK(LWT),INFO(1),IDID,RPAR,IPAR,
  1313.      *   RWORK(LPHI),RWORK(LDELTA),RWORK(LE),
  1314.      *   RWORK(LWM),IWORK(LIWM),
  1315.      *   RWORK(LALPHA),RWORK(LBETA),RWORK(LGAMMA),
  1316.      *   RWORK(LPSI),RWORK(LSIGMA),
  1317.      *   RWORK(LCJ),RWORK(LCJOLD),RWORK(LHOLD),
  1318.      *   RWORK(LS),HMIN,RWORK(LROUND),
  1319.      *   IWORK(LPHASE),IWORK(LJCALC),IWORK(LK),
  1320.      *   IWORK(LKOLD),IWORK(LNS),INFO(10),NTEMP)
  1321. 527   IF(IDID.LT.0)GO TO 600
  1322. C
  1323. C--------------------------------------------------------
  1324. C     THIS BLOCK HANDLES THE CASE OF A SUCCESSFUL RETURN
  1325. C     FROM DDASTP (IDID=1).  TEST FOR STOP CONDITIONS.
  1326. C--------------------------------------------------------
  1327. C
  1328.       IF(INFO(4).NE.0)GO TO 540
  1329.            IF(INFO(3).NE.0)GO TO 530
  1330.              IF((TN-TOUT)*H.LT.0.0D0)GO TO 500
  1331.              CALL DDATRP(TN,TOUT,Y,YPRIME,NEQ,
  1332.      *         IWORK(LKOLD),RWORK(LPHI),RWORK(LPSI))
  1333.              IDID=3
  1334.              T=TOUT
  1335.              GO TO 580
  1336. 530          IF((TN-TOUT)*H.GE.0.0D0)GO TO 535
  1337.              T=TN
  1338.              IDID=1
  1339.              GO TO 580
  1340. 535          CALL DDATRP(TN,TOUT,Y,YPRIME,NEQ,
  1341.      *         IWORK(LKOLD),RWORK(LPHI),RWORK(LPSI))
  1342.              IDID=3
  1343.              T=TOUT
  1344.              GO TO 580
  1345. 540   IF(INFO(3).NE.0)GO TO 550
  1346.       IF((TN-TOUT)*H.LT.0.0D0)GO TO 542
  1347.          CALL DDATRP(TN,TOUT,Y,YPRIME,NEQ,
  1348.      *     IWORK(LKOLD),RWORK(LPHI),RWORK(LPSI))
  1349.          T=TOUT
  1350.          IDID=3
  1351.          GO TO 580
  1352. 542   IF(ABS(TN-TSTOP).LE.100.0D0*UROUND*
  1353.      *   (ABS(TN)+ABS(H)))GO TO 545
  1354.       TNEXT=TN+H
  1355.       IF((TNEXT-TSTOP)*H.LE.0.0D0)GO TO 500
  1356.       H=TSTOP-TN
  1357.       GO TO 500
  1358. 545   CALL DDATRP(TN,TSTOP,Y,YPRIME,NEQ,
  1359.      *  IWORK(LKOLD),RWORK(LPHI),RWORK(LPSI))
  1360.       IDID=2
  1361.       T=TSTOP
  1362.       GO TO 580
  1363. 550   IF((TN-TOUT)*H.GE.0.0D0)GO TO 555
  1364.       IF(ABS(TN-TSTOP).LE.100.0D0*UROUND*(ABS(TN)+ABS(H)))GO TO 552
  1365.       T=TN
  1366.       IDID=1
  1367.       GO TO 580
  1368. 552   CALL DDATRP(TN,TSTOP,Y,YPRIME,NEQ,
  1369.      *  IWORK(LKOLD),RWORK(LPHI),RWORK(LPSI))
  1370.       IDID=2
  1371.       T=TSTOP
  1372.       GO TO 580
  1373. 555   CALL DDATRP(TN,TOUT,Y,YPRIME,NEQ,
  1374.      *   IWORK(LKOLD),RWORK(LPHI),RWORK(LPSI))
  1375.       T=TOUT
  1376.       IDID=3
  1377.       GO TO 580
  1378. C
  1379. C--------------------------------------------------------
  1380. C     ALL SUCCESSFUL RETURNS FROM DDASSL ARE MADE FROM
  1381. C     THIS BLOCK.
  1382. C--------------------------------------------------------
  1383. C
  1384. 580   CONTINUE
  1385.       RWORK(LTN)=TN
  1386.       RWORK(LH)=H
  1387.       RETURN
  1388. C
  1389. C-----------------------------------------------------------------------
  1390. C     THIS BLOCK HANDLES ALL UNSUCCESSFUL
  1391. C     RETURNS OTHER THAN FOR ILLEGAL INPUT.
  1392. C-----------------------------------------------------------------------
  1393. C
  1394. 600   CONTINUE
  1395.       ITEMP=-IDID
  1396.       GO TO (610,620,630,690,690,640,650,660,670,675,
  1397.      *  680,685), ITEMP
  1398. C
  1399. C     THE MAXIMUM NUMBER OF STEPS WAS TAKEN BEFORE
  1400. C     REACHING TOUT
  1401. 610   WRITE (XERN3, '(1P,D15.6)') TN
  1402.       CALL XERMSG ('SLATEC', 'DDASSL',
  1403.      *   'AT CURRENT T = ' // XERN3 // ' 500 STEPS TAKEN ON THIS ' //
  1404.      *   'CALL BEFORE REACHING TOUT', IDID, 1)
  1405.       GO TO 690
  1406. C
  1407. C     TOO MUCH ACCURACY FOR MACHINE PRECISION
  1408. 620   WRITE (XERN3, '(1P,D15.6)') TN
  1409.       CALL XERMSG ('SLATEC', 'DDASSL',
  1410.      *   'AT T = ' // XERN3 // ' TOO MUCH ACCURACY REQUESTED FOR ' //
  1411.      *   'PRECISION OF MACHINE. RTOL AND ATOL WERE INCREASED TO ' //
  1412.      *   'APPROPRIATE VALUES', IDID, 1)
  1413.       GO TO 690
  1414. C
  1415. C     WT(I) .LE. 0.0 FOR SOME I (NOT AT START OF PROBLEM)
  1416. 630   WRITE (XERN3, '(1P,D15.6)') TN
  1417.       CALL XERMSG ('SLATEC', 'DDASSL',
  1418.      *   'AT T = ' // XERN3 // ' SOME ELEMENT OF WT HAS BECOME .LE. ' //
  1419.      *   '0.0', IDID, 1)
  1420.       GO TO 690
  1421. C
  1422. C     ERROR TEST FAILED REPEATEDLY OR WITH H=HMIN
  1423. 640   WRITE (XERN3, '(1P,D15.6)') TN
  1424.       WRITE (XERN4, '(1P,D15.6)') H
  1425.       CALL XERMSG ('SLATEC', 'DDASSL',
  1426.      *   'AT T = ' // XERN3 // ' AND STEPSIZE H = ' // XERN4 //
  1427.      *   ' THE ERROR TEST FAILED REPEATEDLY OR WITH ABS(H)=HMIN',
  1428.      *   IDID, 1)
  1429.       GO TO 690
  1430. C
  1431. C     CORRECTOR CONVERGENCE FAILED REPEATEDLY OR WITH H=HMIN
  1432. 650   WRITE (XERN3, '(1P,D15.6)') TN
  1433.       WRITE (XERN4, '(1P,D15.6)') H
  1434.       CALL XERMSG ('SLATEC', 'DDASSL',
  1435.      *   'AT T = ' // XERN3 // ' AND STEPSIZE H = ' // XERN4 //
  1436.      *   ' THE CORRECTOR FAILED TO CONVERGE REPEATEDLY OR WITH ' //
  1437.      *   'ABS(H)=HMIN', IDID, 1)
  1438.       GO TO 690
  1439. C
  1440. C     THE ITERATION MATRIX IS SINGULAR
  1441. 660   WRITE (XERN3, '(1P,D15.6)') TN
  1442.       WRITE (XERN4, '(1P,D15.6)') H
  1443.       CALL XERMSG ('SLATEC', 'DDASSL',
  1444.      *   'AT T = ' // XERN3 // ' AND STEPSIZE H = ' // XERN4 //
  1445.      *   ' THE ITERATION MATRIX IS SINGULAR', IDID, 1)
  1446.       GO TO 690
  1447. C
  1448. C     CORRECTOR FAILURE PRECEEDED BY ERROR TEST FAILURES.
  1449. 670   WRITE (XERN3, '(1P,D15.6)') TN
  1450.       WRITE (XERN4, '(1P,D15.6)') H
  1451.       CALL XERMSG ('SLATEC', 'DDASSL',
  1452.      *   'AT T = ' // XERN3 // ' AND STEPSIZE H = ' // XERN4 //
  1453.      *   ' THE CORRECTOR COULD NOT CONVERGE.  ALSO, THE ERROR TEST ' //
  1454.      *   'FAILED REPEATEDLY.', IDID, 1)
  1455.       GO TO 690
  1456. C
  1457. C     CORRECTOR FAILURE BECAUSE IRES = -1
  1458. 675   WRITE (XERN3, '(1P,D15.6)') TN
  1459.       WRITE (XERN4, '(1P,D15.6)') H
  1460.       CALL XERMSG ('SLATEC', 'DDASSL',
  1461.      *   'AT T = ' // XERN3 // ' AND STEPSIZE H = ' // XERN4 //
  1462.      *   ' THE CORRECTOR COULD NOT CONVERGE BECAUSE IRES WAS EQUAL ' //
  1463.      *   'TO MINUS ONE', IDID, 1)
  1464.       GO TO 690
  1465. C
  1466. C     FAILURE BECAUSE IRES = -2
  1467. 680   WRITE (XERN3, '(1P,D15.6)') TN
  1468.       WRITE (XERN4, '(1P,D15.6)') H
  1469.       CALL XERMSG ('SLATEC', 'DDASSL',
  1470.      *   'AT T = ' // XERN3 // ' AND STEPSIZE H = ' // XERN4 //
  1471.      *   ' IRES WAS EQUAL TO MINUS TWO', IDID, 1)
  1472.       GO TO 690
  1473. C
  1474. C     FAILED TO COMPUTE INITIAL YPRIME
  1475. 685   WRITE (XERN3, '(1P,D15.6)') TN
  1476.       WRITE (XERN4, '(1P,D15.6)') HO
  1477.       CALL XERMSG ('SLATEC', 'DDASSL',
  1478.      *   'AT T = ' // XERN3 // ' AND STEPSIZE H = ' // XERN4 //
  1479.      *   ' THE INITIAL YPRIME COULD NOT BE COMPUTED', IDID, 1)
  1480.       GO TO 690
  1481. C
  1482. 690   CONTINUE
  1483.       INFO(1)=-1
  1484.       T=TN
  1485.       RWORK(LTN)=TN
  1486.       RWORK(LH)=H
  1487.       RETURN
  1488. C
  1489. C-----------------------------------------------------------------------
  1490. C     THIS BLOCK HANDLES ALL ERROR RETURNS DUE
  1491. C     TO ILLEGAL INPUT, AS DETECTED BEFORE CALLING
  1492. C     DDASTP. FIRST THE ERROR MESSAGE ROUTINE IS
  1493. C     CALLED. IF THIS HAPPENS TWICE IN
  1494. C     SUCCESSION, EXECUTION IS TERMINATED
  1495. C
  1496. C-----------------------------------------------------------------------
  1497. 701   CALL XERMSG ('SLATEC', 'DDASSL',
  1498.      *   'SOME ELEMENT OF INFO VECTOR IS NOT ZERO OR ONE', 1, 1)
  1499.       GO TO 750
  1500. C
  1501. 702   WRITE (XERN1, '(I8)') NEQ
  1502.       CALL XERMSG ('SLATEC', 'DDASSL',
  1503.      *   'NEQ = ' // XERN1 // ' .LE. 0', 2, 1)
  1504.       GO TO 750
  1505. C
  1506. 703   WRITE (XERN1, '(I8)') MXORD
  1507.       CALL XERMSG ('SLATEC', 'DDASSL',
  1508.      *   'MAXORD = ' // XERN1 // ' NOT IN RANGE', 3, 1)
  1509.       GO TO 750
  1510. C
  1511. 704   WRITE (XERN1, '(I8)') LENRW
  1512.       WRITE (XERN2, '(I8)') LRW
  1513.       CALL XERMSG ('SLATEC', 'DDASSL',
  1514.      *   'RWORK LENGTH NEEDED, LENRW = ' // XERN1 //
  1515.      *   ', EXCEEDS LRW = ' // XERN2, 4, 1)
  1516.       GO TO 750
  1517. C
  1518. 705   WRITE (XERN1, '(I8)') LENIW
  1519.       WRITE (XERN2, '(I8)') LIW
  1520.       CALL XERMSG ('SLATEC', 'DDASSL',
  1521.      *   'IWORK LENGTH NEEDED, LENIW = ' // XERN1 //
  1522.      *   ', EXCEEDS LIW = ' // XERN2, 5, 1)
  1523.       GO TO 750
  1524. C
  1525. 706   CALL XERMSG ('SLATEC', 'DDASSL',
  1526.      *   'SOME ELEMENT OF RTOL IS .LT. 0', 6, 1)
  1527.       GO TO 750
  1528. C
  1529. 707   CALL XERMSG ('SLATEC', 'DDASSL',
  1530.      *   'SOME ELEMENT OF ATOL IS .LT. 0', 7, 1)
  1531.       GO TO 750
  1532. C
  1533. 708   CALL XERMSG ('SLATEC', 'DDASSL',
  1534.      *   'ALL ELEMENTS OF RTOL AND ATOL ARE ZERO', 8, 1)
  1535.       GO TO 750
  1536. C
  1537. 709   WRITE (XERN3, '(1P,D15.6)') TSTOP
  1538.       WRITE (XERN4, '(1P,D15.6)') TOUT
  1539.       CALL XERMSG ('SLATEC', 'DDASSL',
  1540.      *   'INFO(4) = 1 AND TSTOP = ' // XERN3 // ' BEHIND TOUT = ' //
  1541.      *   XERN4, 9, 1)
  1542.       GO TO 750
  1543. C
  1544. 710   WRITE (XERN3, '(1P,D15.6)') HMAX
  1545.       CALL XERMSG ('SLATEC', 'DDASSL',
  1546.      *   'HMAX = ' // XERN3 // ' .LT. 0.0', 10, 1)
  1547.       GO TO 750
  1548. C
  1549. 711   WRITE (XERN3, '(1P,D15.6)') TOUT
  1550.       WRITE (XERN4, '(1P,D15.6)') T
  1551.       CALL XERMSG ('SLATEC', 'DDASSL',
  1552.      *   'TOUT = ' // XERN3 // ' BEHIND T = ' // XERN4, 11, 1)
  1553.       GO TO 750
  1554. C
  1555. 712   CALL XERMSG ('SLATEC', 'DDASSL',
  1556.      *   'INFO(8)=1 AND H0=0.0', 12, 1)
  1557.       GO TO 750
  1558. C
  1559. 713   CALL XERMSG ('SLATEC', 'DDASSL',
  1560.      *   'SOME ELEMENT OF WT IS .LE. 0.0', 13, 1)
  1561.       GO TO 750
  1562. C
  1563. 714   WRITE (XERN3, '(1P,D15.6)') TOUT
  1564.       WRITE (XERN4, '(1P,D15.6)') T
  1565.       CALL XERMSG ('SLATEC', 'DDASSL',
  1566.      *   'TOUT = ' // XERN3 // ' TOO CLOSE TO T = ' // XERN4 //
  1567.      *   ' TO START INTEGRATION', 14, 1)
  1568.       GO TO 750
  1569. C
  1570. 715   WRITE (XERN3, '(1P,D15.6)') TSTOP
  1571.       WRITE (XERN4, '(1P,D15.6)') T
  1572.       CALL XERMSG ('SLATEC', 'DDASSL',
  1573.      *   'INFO(4)=1 AND TSTOP = ' // XERN3 // ' BEHIND T = ' // XERN4,
  1574.      *   15, 1)
  1575.       GO TO 750
  1576. C
  1577. 717   WRITE (XERN1, '(I8)') IWORK(LML)
  1578.       CALL XERMSG ('SLATEC', 'DDASSL',
  1579.      *   'ML = ' // XERN1 // ' ILLEGAL.  EITHER .LT. 0 OR .GT. NEQ',
  1580.      *   17, 1)
  1581.       GO TO 750
  1582. C
  1583. 718   WRITE (XERN1, '(I8)') IWORK(LMU)
  1584.       CALL XERMSG ('SLATEC', 'DDASSL',
  1585.      *   'MU = ' // XERN1 // ' ILLEGAL.  EITHER .LT. 0 OR .GT. NEQ',
  1586.      *   18, 1)
  1587.       GO TO 750
  1588. C
  1589. 719   WRITE (XERN3, '(1P,D15.6)') TOUT
  1590.       CALL XERMSG ('SLATEC', 'DDASSL',
  1591.      *  'TOUT = T = ' // XERN3, 19, 1)
  1592.       GO TO 750
  1593. C
  1594. 750   IDID=-33
  1595.       IF(INFO(1).EQ.-1) THEN
  1596.          CALL XERMSG ('SLATEC', 'DDASSL',
  1597.      *      'REPEATED OCCURRENCES OF ILLEGAL INPUT$$' //
  1598.      *      'RUN TERMINATED. APPARENT INFINITE LOOP', -999, 2)
  1599.       ENDIF
  1600. C
  1601.       INFO(1)=-1
  1602.       RETURN
  1603. C-----------END OF SUBROUTINE DDASSL------------------------------------
  1604.       END
  1605.